Wolfram-Powered Polymer Reaction Engineering: A Multiparadigm Approach

Arjun Raghuraman, The Dow Chemical Company

#WolframTechConf

Overview

This talk will cover mathematical modeling of polymerization reactors using the deterministic approach

The polymer of interest are polyether polyols

Graph[{1\[DirectedEdge]2,2\[DirectedEdge]3,4\[DirectedEdge]3,3\[DirectedEdge]5,5\[DirectedEdge]6,6\[DirectedEdge]7},]

Polyether Polyols

Intermediates to polyurethanes; made by polymerization of epoxide monomers, e.g. propylene oxide (PO)

MoleculePlot[Molecule["C1OC1(C)"],ImageSize->Tiny]

Often are star-shaped structures with multiple arms.

An elegant approach to modeling these polymers is to solve for the single-arm distributions, followed by transformation to multi-arm structures.

createMoleculeWithRepeats[n_]:=Module[{repeatUnit,fullStructure},
  repeatUnit="OCC(C)";
  (* Construct the full structure with n repeats *)
  fullStructure=StringJoin["C",StringRepeat[repeatUnit,n],"O[H]"];
  (* Create the molecule *)
  Molecule[fullStructure]
]
repeats=3;
mol1=createMoleculeWithRepeats[repeats];
mol2=MoleculeModify[mol1,{"SetTorsionAngle",Thread[NestList[#+4&,{2,3,4,6},repeats-1]->180]}];
Row[{MoleculePlot[mol2,1],MoleculePlot3D[mol2,1]},"     "]

Integrating Four Paradigms in the Wolfram Language

Symbolic tracking of polymer chains

Leveraging the Paclet repository to generate the population balance (mole balances for all chains and monomer in the reaction system)

Numerical Mathematics: approximate solutions to systems of Ordinary Differential Equations (ODEs) and algebraic equations

Leveraging core language functions: e.g., built-in parametric distributions for initial guesses, and parallel computing capabilities

Symbolic  Tracking of Polymer Chains

The symbolic nature of the Wolfram Language is well suited for application in polymer reaction engineering

can represent the molar concentration of chain with  repeat units

can represent the molar concentration of chain with  repeat units of monomer 1,  repeat units of monomer 2 and end group

In dynamic systems where the concentration of polymer chains varies with time,  can represent the molar concentration of chains of length  as a function of the independent variable

Mechanism of Anionic Ring Opening Polymerization

initiation

propagation

Chain exchange (chain transfer to polymer)

References: (1) Polymer Distribution Models for Polyether Polyols. 2025 North American Symposium for Chemical Reaction Engineering, 17 Feb, Houston TX . (2) Nie, Y., Biegler, L. T., Villa, C. M. & Wassick, J. M. Ind. Eng. Chem. Res. 53, 7434–7446 (2014).

Leveraging the Paclet Repository

General mole balance equation (GMBE):

pts={{0,0},{1,1},{2,-1},{3,1},{4,-2},{2.5,-4},{-0.25,-1}};
Graphics[{Arrow[{{-1,-1},{0.08,-1}}],{Red,BSplineCurve[pts,SplineClosed->True]},Arrow[{{3.55,-1},{3.55+1.08,-1}}],Text["\!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["F", "TI"], 
RowBox[{
StyleBox["j", "TI"], "0"}]], TraditionalForm], "errors" -> {}, "input" -> "F_{j0}", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)",{-1.3,-1}],Text["\!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["F", "TI"], 
StyleBox["j", "TI"]], TraditionalForm], "errors" -> {}, "input" -> "F_j", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)",{3.55+1.3,-1}],Text["\!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["G", "TI"], 
StyleBox["j", "TI"]], TraditionalForm], "errors" -> {}, "input" -> "G_j", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)",{2,-1}]}]

is the molar flow rate of species  at the reactor inlet
 is the molar flow rate of species  at the reactor exit
 is the molar rate of generation of species 
 represents the molar rate of accumulation of species

Reference: Fogler, H. S. Elements of Chemical Reaction Engineering. (Prentice Hall, 2016).

Well-Mixed Batch Reactors

or   where  represents a constant reactor volume and  represents the molar concentration of

Load the compartmental modeling paclet from Robert Nachbar:

PacletInstall["RobertNachbar/CompartmentalModeling"];
Needs["RobertNachbar`CompartmentalModeling`"];

nMax=50(*sets the truncation limit for polyol chains*);rxnNet=KineticCompartmentalModel[{
Table[p[i]+m(Overscript[→, kp])p[i+1],{i,0,nMax-1}],
Table[p[i]+d[j](Overscript[→, ke])p[j]+d[i],{i,0,nMax},{j,0,nMax}]

}//Flatten,t];

rxnNet["Equations"]//Length

Batch Reactor Model & Results

tFinal=4(*hour*);
batchReactorModel[θ1_?NumericQ,θ2_?NumericQ]:=Module[{eqns,ics,vars},
eqns=rxnNet["Equations"];ics=;
vars=rxnNet["Variables"];NDSolve[{eqns,ics}/. {kp->θ1,ke->θ2},vars,{t,0,tFinal}]];

With[{data1=Total/@(TakeDrop[Through[Rest[batchReactorModel[50,(10^6)][[1]]][[All,2]][tFinal]],{1,nMax+1}]\[Transpose]),
data2=Total/@(TakeDrop[Through[Rest[batchReactorModel[50,(10^2)][[1]]][[All,2]][tFinal]],{1,nMax+1}]\[Transpose])},
nlm=NonlinearModelFit[{Range[0,nMax],(data1/(Total[data1]))}\[Transpose],PDF[PoissonDistribution[μ],x],μ,x];
Show[]
]

Well-Mixed Semibatch Reactors

The monomer  is fed into the reactor at a controlled rate for safe operation. Common in industry for exothermic reactions.

Monomer balance:

Polyol chain balance: or  where  is time dependent

The volume of the reaction mixture increases with time.

The compartmental modeling paclet cannot be directly used in this case. Replacement rules were used to modify the equations

Semibatch Reactor Model Results

semiBatchModel[θ1_?NumericQ,θ2_?NumericQ]:=Module[{eqns,ics,unks,feedTime=16,digestTime=2,feedRate},
semiBatchChains=With[
{lhs=Function[D[vol[t] #,t]]/@(Function[Integrate[#,t]]/@rxnNet["Equations"][[Range[2,Length[rxnNet["Equations"]]]]][[All,1]]),
rhs=rxnNet["Equations"][[Range[2,Length[rxnNet["Equations"]]]]][[All,2]]/.x_->vol[t]×x},
Thread[lhs==rhs]
];
feedRate=(81.75/feedTime);pw=Piecewise[{{feedRate,t≤feedTime},
{0,feedTime<t<feedTime+digestTime}}];
;NDSolve[{eqns,ics}/. {kp->θ1,ke->θ2},unks,{t,0,feedTime+digestTime},Method->{"EquationSimplification"->"Residual"}]
]

Copolymers in Batch Reactors

{m1,m2,m3}={Molecule["C1OC1"],Molecule["C1(C)OC1"],Molecule["CO"]};
MapThread[(Labeled[#1,#2]&),{{MoleculePlot3D[m1,ImageSize->Tiny],MoleculePlot3D[m2,ImageSize->Tiny],MoleculePlot3D[m3,ImageSize->Tiny]},{"EO","PO","Methanol"}}]

Let  represent a chain and also the molar concentration of that chain in the reaction mixture where

is the number of PO units independent of the distribution of the PO groups in the chain

is the number of EO units independent of the distribution of the EO groups in the chain

can only take three values: 1 for PO ends, 2 for EO ends and 0 for initiator ends

Copolymers in Batch Reactors

nMax=25;tFinal=20;
rxnNetCopolym=KineticCompartmentalModel[{
p[{0,0},0]+PO(Overscript[→, k10])p[{1,0},1],p[{0,0},0]+EO(Overscript[→, k20])p[{0,1},2],
Table[p[{i,j},1]+EO(Overscript[→, k21])p[{i,j+1},2],{i,nMax},{j,0,nMax}],
Table[p[{i,j},2]+EO(Overscript[→, k22])p[{i,j+1},2],{i,0,nMax},{j,nMax}],
Table[p[{i,j},2]+PO(Overscript[→, k12])p[{i+1,j},1],{i,0,nMax},{j,nMax}],
Table[p[{i,j},1]+PO(Overscript[→, k11])p[{i+1,j},1],{i,nMax},{j,0,nMax}]
}//Flatten,t]//Quiet;//AbsoluteTiming

rxnNetCopolym["Equations"]//Length

modelCopolym[{θ1_?NumericQ,θ2_?NumericQ,θ3_?NumericQ,θ4_?NumericQ,θ5_?NumericQ,θ6_?NumericQ}]:=(modelCopolym[{θ1,θ2,θ3,θ4,θ5,θ6}]=Module[
{eqns,ics,unks},
;
NDSolve[Join[{eqns,ics}]/.{k10->θ1,k20->θ2,k11->θ3,k12->θ4,k21->θ5,k22->θ6},unks,{t,0,tFinal}]
])
modelCopolymResults=modelCopolym[{1,1,1,1,1,1}];
χList=;
ListPlot3D[Flatten[χList,1],AxesLabel->(Style[#,Red]&)/@{"PO", "EO", "χ"},PlotRange->All,ImageSize->350]

Well-Mixed Continuous Stirred Tank Reactors at Steady State

pts={{0,0},{1,1},{2,-1},{3,1},{4,-2},{2.5,-4},{-0.25,-1}};
Graphics[{Arrow[{{-1,-1},{0.08,-1}}],{Red,BSplineCurve[pts,SplineClosed->True]},Arrow[{{3.55,-1},{3.55+1.08,-1}}],Text["\!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["F", "TI"], 
RowBox[{
StyleBox["j", "TI"], "0"}]], TraditionalForm], "errors" -> {}, "input" -> "F_{j0}", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)",{-1.3,-1}],Text["\!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["F", "TI"], 
StyleBox["j", "TI"]], TraditionalForm], "errors" -> {}, "input" -> "F_j", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)",{3.55+1.3,-1}],Text["\!\(\*TemplateBox[<|"boxes" -> FormBox[
SubscriptBox[
StyleBox["G", "TI"], 
StyleBox["j", "TI"]], TraditionalForm], "errors" -> {}, "input" -> "G_j", "state" -> "Boxes"|>,"TeXAssistantTemplate"]\)",{2,-1}]}]

Well-Mixed Continuous Stirred Tank Reactors at Steady State

FindRoot was used to solve this system. Initial guesses were derived from parametric statistical distributions within the core language.

(HyperlinkURL["http://reference.wolfram.com/language/guide/ParametricStatisticalDistributions.html", http://reference.wolfram.com/language/guide/ParametricStatisticalDistributions.html])

cstrModelResults[nMax_]:=With[
,
{eqns=//Flatten},

{igs=Append[
{
Table[n[i],{i,0,nMax}],Table[PDF[PoissonDistribution[(nPOFeed/nH)],x],{x,0,nMax}]
}\[Transpose],
{M,1.}
]
},

{cstrModelResults=FindRoot[eqns,igs]},

ListPlot[{Range[0,nMax],Most[cstrModelResults][[All,2]]}\[Transpose],
]
]

Broader distributions are obtained with a CSTR. The truncation values should be increased.

{cstrModelResults[25],cstrModelResults[100]}//Row

Double Metal Cyanide Catalysis - An Interesting Example

The active catalyst has excess zinc relative to cobalt atoms (~ 2:1)

The cyanocobaltate framework is chemically stable and supports the placement of “extra” zinc within that structure.

The “extra” Zn is the active site.

Empirical Observations on DMC Catalysts

Some unique aspects of DMC catalysis:

The catalyst is more selective to shorter chains than longer chains

The catalyst is inhibited by high concentration of shorter chains

A Langmuir-Hinshelwood mechanism is not sufficient to explain the unusual behavior of these catalysts

Molecular Weight Distributions from DMC Catalysts in CSTRs

Narrow molecular weight distributions are possible even in CSTRs:

In KOH catalysis, K+ can easily find the ends of dormant chains

In the heterogeneous DMC catalyst, the probability of an OH group finding an active Zn center depends on the volume of the random coil shown:

createMoleculeWithRepeats[n_]:=Module[{repeatUnit,fullStructure},
  repeatUnit="OCC(C)";
  (* Construct the full structure with n repeats *)
  fullStructure=StringJoin["C",StringRepeat[repeatUnit,n],"O[Zn]"];
  (* Create the molecule *)
  Molecule[fullStructure]
]
g[molecule_,angle_,repeats_]:=Module[
{torsionChanges,molecules},
torsionChanges=({"SetTorsionAngle",#->angle}&)/@NestList[#+4&,{2,3,4,6},repeats-1];
molecules=Fold[MoleculeModify,molecule,torsionChanges]
]
exampleMol=g[createMoleculeWithRepeats[17],60,17];
backboneIndices=NestList[#+4&,{2,3,4,6},17-1]//Flatten//DeleteDuplicates;
coords=MoleculeValue[exampleMol,{"AtomCoordinates",backboneIndices}];
findGeometricMedian[points_]:=
  NMinimize[Total[Norm[#-{x,y,z}] & /@ points],{x,y,z}][[2,All,2]]
Module[{
scalarCoordinates=QuantityMagnitude[coords],
centralPoint,radius,
g1,g2,f,g
},
centralPoint=findGeometricMedian[scalarCoordinates];
radius=Module[{list=(Norm[#-centralPoint]&)/@scalarCoordinates},
pos=Position[list,Max[list]][[1,1]];
list[[pos]]
];
g1=Graphics3D[{
    {Blue,Tube[scalarCoordinates]},(* Bonds as tubes *)
    {Opacity[0.3],Sphere[centralPoint,radius]} (* Circumscribing sphere *)
  }];
{f,g}={
(MapAt[(Style[#,Red]&),#,-1]&),
(MapAt[(Style[#,Black]&),#,1]&)
};
g2=ListPointPlot3D[Composition[f,g][scalarCoordinates]];
Show[g1,g2]

]

Conclusions

The multi-paradigm and symbolic nature of the Wolfram language offers a powerful platform in polymer reaction engineering.

The compartmental modelling paclet simplifies the process of mathematically defining the polymerization reaction system.

Expansion of this paclet to accommodate well-known reactor configurations  (general mole balance equation) is recommended.

The population balance methodology described in this presentation to solve the reaction system can be limited by the computing power available on a PC. Parallel computing functions alleviate the problem to a certain extent.

Acknowledgements

Carlos Villa (retired from Dow)

The Polyurethanes Process R&D Research team at Dow

Wolfram Support: Robert Nachbar, Louis Breton, Serafine Kovanen

Being trained as a chemist, Wolfram has been instrumental in blurring the interface between chemistry and chemical reaction engineering. Being able to speak both languages is highly beneficial in industrial research.